{ "cells": [ { "cell_type": "markdown", "id": "d8bb6046", "metadata": {}, "source": [ "# Optical Observations of Satellites\n", "\n", "This provides an example of using optical observations of satellites (e.g., from a telescope on the ground) to update the state of the satellite. This is commonly used to update satellite ephemeris, in particular for satellites in MEO and higher orbits that lack easy access to the service volume of the GPS contellation and can be observed for long durations at night without being shadowed by the Earth. The US government and several commercial companies such as [Exoanalytics](https://exoanalytic.com/space-intelligence/) make use of these techniques to maintain the ephemeris on these satellites.\n", "\n", "This example uses a batch least-squares technique for the fitting. This can also be performed with sequential filters such as an Extended or Unscented Kalman Filter.\n", "\n", "This example updates the initial state vector to mimimizes the difference between observed and predicted line of site to the satellite for all observations:\n", "\n", "$$\n", "\\hat{\\mathbf x}_0\n", "=\n", "\\arg\\min_{\\mathbf x_0\\in\\mathbb R^6}\n", "\\sum_{k=1}^{N}\n", "\\left\\|\n", "E_k\\Big(\\hat{\\mathbf u}_m(\\alpha_k,\\delta_k)-\\hat{\\mathbf u}(\\phi(t_k,t_0;\\mathbf x_0),\\mathbf r_{0}(t_k))\\Big)\n", "\\right\\|_{R^{-1}}^{2}\n", "$$\n", "\n", "\n", "where \n", "$$\n", "\\hat{\\mathbf u}_m(\\alpha,\\delta)=\n", "\\begin{bmatrix}\n", "\\cos\\delta\\cos\\alpha\\\\\n", "\\cos\\delta\\sin\\alpha\\\\\n", "\\sin\\delta\n", "\\end{bmatrix}\n", "$$\n", "is the measurement ($\\alpha$ is right ascension and $\\delta$ is declination),\n", "\n", "$$\n", "\\hat{\\mathbf u}(\\mathbf x_k,\\mathbf r_0)=\n", "\\frac{\\mathbf r_k-\\mathbf r_0}{\\|\\mathbf r_k-\\mathbf r_0\\|}\n", "$$\n", "is the predicted unit vector pointing from the observer at $\\mathbf r_0$ to the satellite at $\\mathbf r_k$ \n", "$$\n", "\\mathbf x_k=\\phi(t_k,t_0;\\mathbf x_0)\n", "$$\n", "is the propagated state at time $t_k$, given initial state $\\mathbf x_0$\n", "\n", "$E$ is just a rotation into the local tangent plane for each observation (to avoid singularities and keep measurement errors constant in local measurement plane)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "a53d5d64", "metadata": {}, "outputs": [], "source": [ "import satkit as sk\n", "import numpy as np\n", "import math as m\n", "import datetime\n", "\n", "# Take our observer to be the Space Surveillance Telescope (SST)\n", "# in Exmouth, Western Australia\n", "sst_lat = -22.675\n", "sst_lon = 114.235\n", "sst_alt = 0.5 # km\n", "sst = sk.itrfcoord(latitude_deg=sst_lat, longitude_deg=sst_lon, alt_m=sst_alt*1e3)\n", "\n", "# Lets look at INTELSAT 19 (IS-19), a geostationary communications satellite.\n", "# It is currently located at 166E longitude, so it should be visible from the SST.\n", "# We can get the TLE for INTELSAT 19 from Celestrak:\n", "# https://celestrak.com/NORAD/elements/geo.txt\n", "tle = sk.TLE.from_lines([\n", " \"INTELSAT 19 (IS-19)\",\n", " \"1 38356U 12030A 26045.51850411 -.00000051 00000+0 00000+0 0 9998\",\n", " \"2 38356 0.0150 97.8441 0002877 233.5096 165.8257 1.00272325 49221\"\n", "])\n", "epoch = tle.epoch\n", "\n", "# Lets make the initial \"true\" state the tle-predicted state at epoch\n", "pteme, vteme = sk.sgp4(tle, epoch)\n", "q = sk.frametransform.qteme2gcrf(epoch)\n", "pgcrf = q * pteme\n", "vgcrf = q * vteme\n", "state0 = np.hstack((pgcrf, vgcrf))\n", "print(epoch)\n", "begin_time = epoch\n", "end_time = epoch + sk.duration(days=2)\n", "settings = sk.propsettings()\n", "settings.precompute_terms(begin_time, end_time)\n", "# True state over 1 day\n", "truth = sk.propagate(state0, begin_time, end_time, propsettings=settings)\n", "\n", "# Assume position knowledge to 10 km, velocity to 5 cm/s\n", "pos_noise = 10e3\n", "vel_noise = 0.05\n", "\n", "# initial state estimate\n", "state0_est = state0 + np.random.normal(scale=[pos_noise]*3 + [vel_noise]*3)\n", "state0_est_prior = state0_est.copy()\n", "\n", "# sample times for measurements. Every 5 minues for 6 hours beginning at local 9pm and going to local 4am\n", "sample_start = sk.time(2026, 2, 15, 20, 0, 0) # local 9pm at SST\n", "sample_end = sample_start + sk.duration(hours=7)\n", "sample_times = [sample_start + sk.duration(minutes=5)*i for i in range(7*60//5)]\n", "\n", "# Truth at first sample time\n", "state0_truth = truth.interp(sample_times[0])\n", "\n", "# the optical measurement is a line-of-sight vector from observer to satellite in the inertial frame\n", "sst_pos = [sk.frametransform.qitrf2gcrf(t) * sst.vector for t in sample_times]\n", "sat_state = [truth.interp(t) for t in sample_times]\n", "los_meas = [sat_state[idx][0:3] - sst_pos[idx] for idx in range(len(sample_times))]\n", "# Normalize the line-of-sight measurements to get unit vectors\n", "los_meas = [los / np.linalg.norm(los) for los in los_meas]\n", "\n", "# add noise to the measurements. Assume 30 microradians of angular noise\n", "# Note: I made up 30 microradians (6 arcseconds).\n", "ang_noise = 30e-6\n", "\n", "def add_ang_noise(u, ang_noise):\n", " # add noise in a random direction perpendicular to the line of sight\n", " # first get a random vector\n", " rand_vec = np.random.normal(size=3)\n", " # make it perpendicular to the line of sight\n", " rand_vec -= rand_vec.dot(u) * u\n", " rand_vec /= np.linalg.norm(rand_vec)\n", " # rotate the line of sight by the noise angle in the direction of the random vector\n", " theta = np.random.normal(scale=ang_noise)\n", " u_noisy = u * m.cos(theta) + np.cross(rand_vec, u) * \\\n", " m.sin(theta) + rand_vec * rand_vec.dot(u) * (1 - m.cos(theta))\n", " return u_noisy\n", "\n", "los_meas_noisy = [add_ang_noise(los, ang_noise) for los in los_meas]\n", "\n", "\n", "# Now, make up an initial state estimate that we are going to refine with measurements\n", "state0_est = truth.interp(sample_times[0]) + np.random.normal(scale=[pos_noise]*3 + [vel_noise]*3)\n", "\n", "# Our initial state error .. record for later comparison\n", "state0_est_prior = state0_est.copy()\n", "\n", "def unit_vector_jacobian(rho):\n", " \"\"\"\n", " Jacobian of unit vector with respect to the original vector\n", " \"\"\"\n", " rho_norm = np.linalg.norm(rho)\n", " u = rho / rho_norm\n", " return (np.eye(3) - np.outer(u, u)) / rho_norm\n", "\n", "\n", "def tangent_basis(u):\n", " \"\"\"\n", " RA/Dec tangent basis at predicted direction\n", " \"\"\"\n", " ux, uy, uz = u\n", " alpha = np.arctan2(uy, ux)\n", " delta = np.arcsin(uz)\n", "\n", " e_alpha = np.array([-np.sin(alpha), np.cos(alpha), 0.0])\n", " e_delta = np.array([\n", " -np.sin(delta)*np.cos(alpha),\n", " -np.sin(delta)*np.sin(alpha),\n", " np.cos(delta)\n", " ])\n", " return np.vstack((e_alpha, e_delta))\n", "\n", "# OK, lets do it!\n", "# 3 iterations of Gauss-Newton least squares\n", "for n in range(0,3):\n", " print(f\"Iteration {n}\")\n", "\n", " # Propagate state to get state transition matrix at each time\n", " res = sk.propagate(state0_est, sample_times[0], sample_times[-1], output_phi=True, propsettings=settings)\n", " [state,phi] = zip(*[res.interp(t, output_phi=True) for t in sample_times])\n", "\n", " residuals = []\n", " J_list = []\n", "\n", " for idx, t in enumerate(sample_times):\n", " # Get the predicted line of sight measurement at this time\n", " p_sat = state[idx][0:3]\n", " p_obs = sst_pos[idx]\n", " los_vec = p_sat - p_obs\n", " los_pred = los_vec / np.linalg.norm(los_vec)\n", "\n", " # Get the tangent basis in predicted direction\n", " E = tangent_basis(los_pred)\n", "\n", " # Get residuals\n", " r = E @ (los_meas_noisy[idx] - los_pred)\n", " residuals.append(r)\n", "\n", " # IMPORTANT: Jacobian must be evaluated on un-normalized LOS vector\n", " U = unit_vector_jacobian(los_vec)\n", " # A is our state transition matrix for this time,\n", " # which maps changes in initial state to changes in state at this time.\n", " # We only care about the position part of the state\n", " # since the measurement is only a function of position.\n", " A = phi[idx][0:3, :]\n", " Hk = - E @ U @ A\n", " J_list.append(Hk)\n", "\n", " # OK, now do the least squares solve to get state update\n", " r = np.hstack(residuals)\n", " H = np.vstack(J_list)\n", " # Solve for state update using least squares\n", " dx = np.linalg.lstsq(H, r, rcond=None)[0]\n", " # Adjust the initial state estimate by the state update\n", " state0_est -= dx\n", "\n", "print(f\"Initial State Error: {np.linalg.norm(state0_est_prior - state0_truth)/1e3:.3f} km, {np.linalg.norm(state0_est_prior[3:6] - state0_truth[3:6]):.3f} m/s\")\n", "print(f'Final State Error: {np.linalg.norm(state0_est - state0_truth)/1e3:.3f} km, {np.linalg.norm(state0_est[3:6] - state0_truth[3:6]):.3f} m/s')" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.14.0" } }, "nbformat": 4, "nbformat_minor": 5 }